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^^ Abstract. The focus of this paper is on the analysis of the Conjugate 

^.^ Gradient method applied to a non-symmetric system of linear equations, 

.^fH arising from a Fast Fourier Transform-based homogenization method 

I "^ due to Moulinec and Suquet [T] . Convergence of the method is proven 

^H by exploiting a certain projection operator reflecting physics of the un- 

ri derlying problem. These results are supported by a numerical example, 

demonstrating significant improvement of the Conjugate Gradient-based 
scheme over the original Moulinec-Suquet algorithm. 
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^. The last decade has witnessed a rapid development in advanced experimen- 

r^ tal techniques and niodehng tools for microstructural characterization, typically 

provided in the form of pixel- or voxel-based geometry. Such data now allow 
for the design of bottom-up predictive models of the overall behavior for a wide 
range of engineering materials. Of course, such step necessitates the develop- 
ment of specialized algorithms, capable of handling large-scale voxel-based data 
in an efficient manner. In the engineering community, perhaps the most success- 
5^ ful solver meeting these criteria was proposed by Moulinec and Suquet in [T]. 

The algorithm is based on the Neumann series expansion of the inverse of an 
operator arising in the associated Lippmann-Schwinger equation and exploits 
the Fast Fourier Transform to evaluate the action of the operator efficiently for 
voxel-based data. In our recent work [2], we have offered a new approach to the 
Moulinec-Suquet scheme, by exploiting the trigonometric collocation method due 
to Saranen and Vainikko [3j. Here, the Lippman-Schwinger equation is projected 
to a space of trigonometric polynomials to yield a non-symmetric system of lin- 
ear equations, see Section [2] below. Quite surprisingly, numerical experiments 
revealed that the system can be efficiently solved using the standard Conjugate 



Gradient algorithni. The analysis of this phenomenon, as presented in Section [3] 
is at the heart of this contribution. The obtained resuhs are further supported 
by a numerical example in Section HI and summarized in Section [Sj 

The following notation is used throughout the paper. Symbols a, a and A 
denote scalar, vector and second-order tensor quantities, respectively, with Greek 
subscripts used when referring to the corresponding components, e.g. AajS. The 
outer product of two vectors is denoted as a® a, whereas a-b oi A-b represents 
the single contraction between vectors (or tensors). A multi-index notation is 
employed, in which K^ with N = {Ni, ...,Nd) represents rN^x-xn^ ^j^j |jy| 
abbreviates Oa^i -^a- Block matrices are denoted by capital letters typeset in a 
bold serif font, e.g. A e ^dxdxNxN ^ g^^^^j ^-^^ superscript and subscript indexes 

are used to refer to the components, such that A = [A^a*]^'^]^ ^ with 



N _ ] < ^ rj,d 



N N 

\ 2 " 2 ' 

Sub-matrices of A are denoted as 

al3=[/^a0\ eM ,A = LAa/3 J„_^^i_..._<ielR 

for a, /3 = 1, . . . ,d and k,m e Z^ . Analogously, the block vectors are denoted 
by lower case letters, e.g. e e M.'^^-'^ and the matrix-by- vector multiplication is 
defined as 

[Ae]^=S S A^^e^eM'^x^, (1) 

with a = 1, . . . , d and k e Z^ . 

2 Problem setting 

Consider a composite material represented by a periodic unit cell 



y= Yl{-Ya,Y^) 



nd 
— -^ on J- a ) ^ •^'>- ■ 



In the context of linear electrostatics, the associated unit cell problem reads as 

V X e(a;) = 0, V-j(a;) = 0, j{x) = L{x) ■ e{x), xey (2) 

where e is a J^-periodic vectorial electric field, j denotes the corresponding vector 
of electric current and JD is a second-order positive-definite tensor of electric 
conductivity. In addition, the field e is subject to a constraint 

(eix)):=^^j eix)dx = e^, (3) 



where |3^| denotes the d-dimensional measure of y and e" # a prescribed 
macroscopic electric field. 

The original problem ([2])-(l3| is then equivalent to the periodic Lippmann- 
Schwinger integral equation, formally written as 

e{x) +jr{x- y; i°) • (L{y) - L°) ■ e{y) dy = e°, xe y, (4) 

where L'^ e M'*^'* denotes a homogeneous reference medium. The operator 
r{x, i") is derived from the Green's function of the problem (|2])-(l3| with 
L{x) = L and can be simply expressed in the Fourier space 

Operator / = /(fe) stands for the Fourier coefficient of f{x) for the fc-th fre- 
quency given by 



f{k) = f{x)ip_k{x)dx, (pk{x) = \y\ 2 exp ivr V -^^ , 



(6) 



"i" is the imaginary unit (i^ = —1). We refer to |2|4| for additional details. 
Note that the linear electrostatics serves here as a model problem; the frame- 
work can be directly extended to e.g. elasticity [5], (visco-)plasticity |;6 or to 
multiferroics [7]. 

2.1 Discretization via trigonometric collocation 

The numerical solution of the Lippmann-Schwinger equation is based on a dis- 
cretization of a unit cell y into a regular periodic grid with A^i x • • • x Nd nodal 
points and grid spacings h = {2Yi/Ni, . . . , 2Yd/Nd)- The searched field e in Q 
is approximated by a trigonometric polynomial e''^ in the form (cf. [Sj Chapter 
10]) 

e{x) « e'^ix) = 2 eVfc(a;), x e y, (7) 

where e = (e^)^^!,...,^ designates the Fourier coefficients defined in ^. No- 
tice that the trigonometrical polynomials are uniquely determined by a regular 
grid data, which makes them well-suited to problems with pixel- or voxel-based 
computations. 

The trigonometric collocation method is based on the projection of the 
Lippmann-Schwinger equation Q onto the space of the trigonometric polyno- 
mials 

r^ = { 2 Cfc(^fc,cfc6c}, (8) 



leading to a to linear system in the form, cf. [2] 

(l + B)e = e«, B = F-ifF(L-L''), (9) 



femsZ™ 



where e = (eS)^^^_ _^ e M''^-^ is the unknown vector, I = [Sai3Skm]a,,f3^i,...,d 
^dxdxNxN jg ^YiQ identity matrix, expressed as the product of the Kronecker 
delta functions So./} and 4m, and e" = (e^)^fi'^ ^ e M''^^. 

All the matrices in ^ exhibit a block-diagonal structure. In particular. 



[4m r^^ J 



k-.meZ _ ^^ , fc^Tfc,meZ~ .n re- .n ik,mel 

a,l3^1,...,d ' 



re I fcm,-|'=^'"SZ [x |0 -]'='■ 

L4mL„/3 Ja^^ = l,...,d' ■- - L'^*""'- «/3Ja,/3 = l„ 



with f^^ = fa,p{k;L"), L^^ = L„^(fc) and (L0)„^ = L%. The matrix F imple- 
mcnts the Discrete Fourier Transform and is defined as 



,fc,msZ~ ^km _ |y|- 






with F^ representing the inverse transform. 

It follows from Eq. Q that the cost of multiplication by B is dominated by 
the action of F and F^ , which can be performed in 0(|Ar| log |7V|) operations 
by the Fast Fourier Transform techniques. This makes the system ^ well-suited 
for applying some iterative solution technique. In particular, the original Fast 
Fourier Transform-based Homogenization scheme formulated by Moulinec and 
Suquet in [T] is based on the Neumann expansion of the matrix inverse (I -I- B)^^, 
so as to yield the rn-th iterate in the form 



e 



(m) 



E(-Bre''. (11) 



j=0 



As indicated earlier, our numerical experiments [2] suggest that the system can 
be efficiently solved using the Conjugate Gradient method, despite the non- 
symmetry of B evident from (|9]). This observation is studied in more detail in 
the next Section. 



3 Solution by the Conjugate Gradient method 

We start our analysis with recasting the system ([9| into a more convenient form, 
by employing a certain operator and the associated sub-space introduced later. 
Note that for simplicity, the reference conductivity is taken as L" = Al. 

Definition 1. Given A > 0, we define operator Pg = AF^ FF and associated 
sub-space as 

f ={P£xforxeM'*^^}cM'^^^. 

Lemma 1. The operator P^ is an orthogonal projection. 



Proof. First, we will prove that P^ is projection, i.e. P^ = P^. Since F is a 
unitary matrix, it is easy to see that 

p2 = (AF-ifF)(AF-ifF) = F-i(Af)2F. (12) 

Hence, in view of the block-diagonal character of f, it it sufficient to prove the 
projection property of sub-matrices (AT)'"'" only. This follows using a simple 
algebra, recall Eq. ([s]): 

The orthogonality of Pg now follows from 

P* = (^AF-ifF)* = AF*r* (F-i)* = AF-ifF = P^, 

according to a well-known result of linear algebra, e.g. Proposition 1.8 in [8^. □ 

Remark 1. It follows from the previous results that the subspace £ collects the 
non-zero coefficients of trigonometric polynomials 7"''^ with zero rotation, which 
represent admissible solutions to the unit cell problem defined by (pi) . Note that 
the orthogonal space £ contains the trigonometric representation of constant 
fields, cf. H Section 12.7]. 

Lemma 2. The solution e to the linear system ([9]) admits the decomposition 
e = e + &£ , with eg e £ satisfying 

P£Le£ + P£Le'' = 0. (13) 

Proof. As e 6 R'^^^, Lemma [l] ensures that it can be decomposed into two 
orthogonal parts e^ = P^e and e^x = (' — Ps)^- Substituting this expression 
into Q, and using the identity B = AF^^FF (-^ — I), we arrive at 

-PeLes + e£± + ^PfLe^x = e". (14) 

Since e° e £-*", we have e^x = e" and the proof is complete. □ 

With these auxiliary results in hand, we are in the position to present our 
main result. 

Proposition 1. The non-symmetric system of linear equations Q is solvable 
by the Conjugate Gradient method for an initial vector etQ\ = e" -I- e with e e £. 
Moreover, the sequence of iterates is independent of the parameter A. 

Proof (outline). It follows from Lemma [2] that the solution to dol admits yet 
another, optimization-based, characterization in the form 



e -I- arg min 

eeS 



-(Le,e)^,,„ + (Le", 



(15) 



The residual corresponding to the initial vector ejo) equals to 

1_ . n 1, 



r(o) 



e" - (I + B) (e° + e) = --P^Le" - ^PfLe e £ 



A^ A 



It can be verified that the subspace £ is B-invariant, thus (l + B)f a £. Therefore, 
the Krylov subspace 

jr™(l + B, r(o)) = span{r(o), (I + B)r(o), . . . , (I + B)'"r(o)} c £ 

for arbitrary m e N. This implies that the residual ri^\ and the Conjugate Gra- 
dient search direction p,,„-) at the m-th iteration satisfy ri^\ e £ and P(„) e £. 
Since B is symmetric and positive-definite on £", the convergence of CG algo- 
rithm now follows from standard arguments, e.g. Theorem 6.6 in 0. Observe 
that different choices of A generate identical Krylov subspaces, thus the sequence 
of iterates is independent of A. □ 

Remark 2. Note that it is possible to show, using direct calculations based on the 
projection properties of Pg, that the Biconjugate Gradient algorithm produces 
exactly the same sequence of vectors as the Conjugate Gradient method, see [S]. 



4 Numerical example 

To support our theoretical results, we consider a three-dimensional model prob- 
lem of electric conduction in a cubic periodic unit cell y = Y[a=i{~\^ \)^ '^^P' 
resenting a two-phase medium with spherical inclusions of 25% volume fraction. 
The conductivity parameters are defined as 



L{x) 



/ 1 0.2 0.2\ 
I 0.2 1 0.2 I 
^ \0.2 0.2 1 / 



otherwise 



where p > denotes the contrast of phase conductivities. We consider the 
macroscopic field e° = [1,0,0] and discretize the unit cell with TV = [n,ri,n] 
node^ The conductivity of the homogeneous reference medium L'^ e M'^^'' is 
parametrized as 

L° = \I, \ = l-uj + puj, (16) 

where w « 0.5 delivers the optimal convergence of the original Moulinec-Suquet 
Fast-Fourier Transform-based Homogcnization (FFTH) algorithm [T]. 

We first investigate the sensitivity of Conjugate Gradient (CG) algorithm to 



the choice of reference medium. The results appear in Fig. 1(a) plotting the 



In particular, n was taken consequently as 16, 32, 64, 128 and 160 leading up to 
3 • 160^ = 12.2 X 10*^ unknowns 



relative number of iterations for CG against the conductivity of the reference 
medium parametrized by w, recall Eq. (16). As expected, CG solver achieve a 



significant improvement over FFTH method as it requires about 40% iterations 
of FFTH for a mildly-contrasted composite down to 4% for g = 10'^. The minor 
differences visible especially for p = 10"^ can be therefore attributed to accumula- 
tion of round-off errors. These observations fully confirm our theoretical results 
presented earlier in Section [3) 
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Fig. 1. (a) Relative number of iterations as a function of the reference medium param- 
eter uj and (b) computational time as a function of the number of unknowns. 



In Fig. 1(b) we present the total computational tima^las a function of the 
number of degrees of freedom and the phase ratio p. The results confirm that 
the computational times scales linearly with the increasing number of degrees of 
freedom for both schemes for a fixed p [5] . The ratio of the computational time 
for GG and FFTH algorithms remains almost constant, which indicates that the 
cost of a single iteration of GG and FFTH method is comparable. 

In addition, the memory requirements of both schemes are also comparable. 
This aspect represents the major advantage of the short-recurrence CG-based 
scheme over alternative schemes for non-symmetric systems, such as GMRES. 
Finally note that finer discretizations can be treated by a straightforward parallel 
implementation. 



5 Conclusions 



In this work, we have proven the convergence of Conjugate Gradient method for 
a non-symmetric system of linear equations arising from periodic unit cell ho- 

'^ The problem was solved with a Matlab® in-house code on a machine Intel® 
Core'^^2 Duo 3 GHz CPU, 3.28 GB computing memory with Debian hnux 5.0 
operating system. 



mogenization problem and confirmed it by numerical experiment. The important 
conclusions to be pointed out are as follows: 

1. The success of the Conjugate Gradient method follows from the projection 
properties of operator P^ introduced in Definition [T] which reflect the struc- 
ture of the underlying physical problem. 

2. Contrary to all available extensions of the FFTH scheme, the performance 
of the Conjugate Gradient-based method is independent of the choice of 
reference medium. This offers an important starting point for further im- 
provements of the method. 

Apart from the already mentioned parallelization, performance of the scheme 
can further be improved by a suitable preconditioning procedure. This topic is 
currently under investigation. 
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